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Abstract 

Understanding the dynamic stability of bodies in frictional contact steadily sliding one over the 
other is of basic interest in various disciplines such as physics, solid mechanics, materials science 
and geophysics. Here we report on a two-dimensional linear stability analysis of a deformable 
solid of a finite height H, steadily sliding on top of a rigid solid within a generic rate-and-state 
friction type constitutive framework, fully accounting for elastodynamic effects. We derive the 
linear stability spectrum, quantifying the interplay between stabilization related to the frictional 
constitutive law and destabilization related both to the elastodynamic bi-material coupling between 
normal stress variations and interfacial slip, and to finite size effects. The stabilizing effects related 
to the frictional constitutive law include velocity-strengthening friction (i.e. an increase in frictional 
resistance with increasing slip velocity, both instantaneous and under steady-state conditions) and 
a regularized response to normal stress variations. We first consider the small wave-number k limit 
and demonstrate that homogeneous sliding in this case is universally unstable, independently of 
the details of the friction law. This universal instability is mediated by propagating waveguide¬ 
like modes, whose fastest growing mode is characterized by a wave-number satisfying kH 0(1) 
and by a growth rate that scales with H We then consider the limit kH —?• oo and derive 
the stability phase diagram in this case. We show that the dominant instability mode travels at 
nearly the dilatational wave-speed in the opposite direction to the sliding direction. In a certain 
parameter range this instability is manifested through unstable modes at all wave-numbers, yet 
the frictional response is shown to be mathematically well-posed. Instability modes which travel at 
nearly the shear wave-speed in the sliding direction also exist in some range of physical parameters. 
Previous results obtained in the quasi-static regime appear relevant only within a narrow region 
of the parameter space. Finally, we show that a finite-time regularized response to normal stress 
variations, within the framework of generalized rate-and-state friction models, tends to promote 
stability. The relevance of our results to the rupture of bi-material interfaces is briefly discussed. 
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1. Background and motivation 


The dynamic stability of steady-state homogeneous sliding between two macroscopic bodies in 
frictional contact is a basic problem of interest in various scientific disciplines such as physics, solid 
mechanics, materials science and geophysics. The emergence of instabilities may give rise to rich 


dynamics and play a dominant role in a broad range of frictional phenomena ( 

Ruina 

T983| Ben- 

Zion, 2001 iScholz, 2002 

Ben-Zion 

2008 

Gerde and Marder 

2001; Ibrahim, 1994a|b[ Di Bartolomeo 
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et al., 2010; Tonazzi et al., 2013 Baillet et al., 2005; Behrendt et al., 2011 Meziane et al. 2007). The 


response of a frictional system to spatiotemporal perturbations, and the accompanying instabilities, 
are governed by several physical properties and processes. Generally speaking, one can roughly 
distinguish between bulk effects (i.e. the constitutive behavior and properties of the bodies of 
interest, their geometry and the external loadings applied to them) and interfacial effects related 
to the frictional constitutive behavior. The ultimate goal of a theory in this respect is to identify the 
relevant physical processes at play, to quantify the interplay between them through properly defined 
dimensionless parameters and to derive the stability phase diagram in terms of these parameters, 
together with the growth rate of various unstable modes. 

As a background and motivation for what will follow, we would like first to briefly discuss the 
various players affecting the stability of frictional sliding, along with stating some relevant results 
available in the literature. Focusing first on bulk effects, it has been recognized that when consid¬ 
ering isotropic linear elastic bodies, there is a qualitative difference between sliding along interfaces 
separating bodies made of identical materials and interfaces separating dissimilar materials. In the 
former case, there is no coupling between interfacial slip and normal stress variations, while in 


the latter case — due to broken symmetry 

— such coupling exists (Comninou 

1977a 

b[ Comni- 

nou and Schmueser, 1979; Weertman 1980 

Andrews and Ben-Zion , |1997[ Ben-Zi 

on anc 

Andrews 

1998[ Adams 

2000t Cochard and Rice 2000 Rice et al. 2001[ Ranjith and Rice 

2001 

Cerde and| 

Harder, 2001 

Adda-Bedia and Ben Amar 

2003t Ampuero and Ben-Zion, 2008 

). This coupling 


may lead to a reduction in the normal stress at the interface and consequently to a reduction in 
the frictional resistance. Hence, bulk material contrast (i.e. the existence of a bi-material interface) 
potentially plays an important destabilizing role in the stability of frictional sliding. Another class 
of bulk effects is related to the hnite geometry of any realistic sliding bodies and the type of loading 
applied to them (e.g. velocity or stress boundary conditions). To the best of our knowledge, the 


latter effects are significantly less explored in the literature (but see Rice and Ruina (1983); Ranjith 


(2009 2014)). 


In relation to interfacial effects, it has been established that sliding along a bi-material interface 
described by the classical Coulomb friction law, T = af (r is the local friction stress, a is the local 
compressive normal stress and / is a constant friction coefficient), is unstable against perturbations 
at all wavelengths and irrespective of the value of the friction coefficient, when the bi-material 
contrast is such that the generalized Rayleigh wave exists (Ranjith and Rice, 2001). The latter 


is an interfacial wave that propagates along frictionless bi-material interfaces, constrained not to 
feature opening (Weertman, 1963 Achenbach and Epstein, jl967 Adams, 1998 Ranjith and Rice 


2001). It is termed the generalized Rayleigh wave because it coincides with the ordinary Rayleigh 


wave when the materials are identical and it exists when the bi-material contrast is not too large. In 


fact, the response to perturbations in this case is mathematically ill-posed ( 

Renardy 

19921 

Adams 

1995 

Martins and Simoes, 1995t Martins et al. 1995;|Sim6es and Martins, 1998[ Ranjith and Rice 

2001 

. Ill-posedness, which is a stronger condition than instability (i.e. all perturbation modes can 


be unstable, yet a problem can be mathematically well-posed), will be discussed later. It has been 
then shown that replacing Coulomb friction by a friction law in which the friction stress r does not 
respond instantaneously to variations in the normal stress a, but rather approaches r = af over 


a finite time scale, can regularize the problem, making it mathematically well-posed (Ranjith and 


Rice 2001). 


Subsequently, the problem has been addressed within the constitutive framework of rate-and- 
state friction models, where the friction stress depends both on the slip velocity and the structural 
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state of the interface. Within this framework (Dieterich 

1978 

19791 

Ruina 

1983 

Rice and Ruina 

1983[ Heslot et ah, 1994| Maronel 1998; Berthoud et a 

.1 1999[ Baumberger anc 

Berthoudl |1999 


Baumberger and Caroli, 2006), the simplest version of the friction stress takes the form T = af{(j), v), 


where v is the difference between the local interfacial slip velocities of the two sliding bodies and (p 
is a dynamic coarse-grained state variabl^ Under steady-state sliding conditions the state variable 
(j) attains a steady-state value determined by v, po^v). Within such a constitutive framework, the 
most relevant physical quantities for the question of stability, which will be extensively discussed 
below, are the instantaneous response to variations in the slip velocity, dj{(j),v) (the so-called 
“direct effect”), and the variation of the steady-state frictional strength with the slip velocity, 

Note that here and below we use the following shorthand notation: 


2001 ) 


dyf {4>o{v), v) (Rice et al 
dv^ii and dv = -§^. 

Previous studies have argued that an instantaneous strengthening response, i.e. the experi¬ 
mentally well-established positive direct effect dvf{4>, v) > 0 (which is associated with thermally 
activated rheology (Baumberger and Caroli 2006)), is sufficient to give rise to the existence of 


a quasi-static range of response to perturbations at sufficiently small slip velocities (Rice et al. 


2001). The existence of such a quasi-static regime is non-trivial (e.g. it does not exist for Coulomb 


friction); it implies that when very small slip velocities are of interest, one can reliably address the 
stability problem in the framework of quasi-static elasticity (i.e. omitting inertial terms to begin 
with), rather than considering the full — and more difficult — elastodynamic problem and then take 
the quasi-static limit. Within such a quasi-static framework, it has been shown that dyf((j),v) >0 
can lead to stable response against sufficiently short wavelength perturbations, even if the inter¬ 
face is velocity-weakening in steady-state, dvf{4>o{v),v) < 0 (Rice et al., 2001). Furthermore, it 


has been shown that sufficiently strong velocity-strengthening, dyf{4>Q{v),v)>0, can overcome the 
destabilizing bi-material effect, leading to the stability of perturbations at all wavelengths in the 
quasi-static limit ( Rice et akf 2001). 

Despite this progress, several important questions remain open. First, to the best of our 
knowledge the fully elastodynamic stability analysis of bi-material interfaces in the framework of 
rate-and-state friction models has not been performed. This is important since the quasi-static 
regime — when it exists — is expected to be valid only for very small slip velocities (as was argued 


m 


Rice et al. (2001) and will be explicitly shown below). Second, a very recent compilation of a 


large set of experimental data for a broad range of materials (]Bar-Sinai et al. 2014) has revealed 


that dry frictional interfaces generically become velocity-strengthening over some range of slip 


velocities 

Weeksl 

1993 

Marone and Scholz, 1988 

Marone et al. 

1991 

Kato 

2003 

Shibazaki and 

lio 

2003 

Bar Sinai et al. 

2012 

Hawthorne and Rubin 

2013 

Bar-Sinai et al. 

2013 

2015 

). In 


other cases, frictional interfaces are intrinsically velocity-strengthening (Perfettini and Ampuero 


2008, Noda and Shimamoto 


2009 

Ikari et al. 

2009 

2013 


As velocity-strengthening friction 
is expected to play a stabilizing role in the stability of frictional sliding, there emerges a basic 
question about the interplay between the stabilizing velocity-strengthening friction effect and the 
destabilizing bi-material effect, when elastodynamics is fully taken into account. Finally, in almost 
all of the previous studies we are aware of, the sliding bodies were assumed to be infinite (but see, 


for example. Rice and Ruina (1983); Ranjith (2009, 2014)). Yet, realistic sliding bodies are of finite 
extent and the interaction with the boundaries may be of importance. 

To address these issues we analyze in this paper the linear stability of a deformable solid of 
height H steadily sliding on top of a rigid solid within a generic rate-and-state friction constitutive 


^In principle there can be more than one internal state variables, but we do not consider this possibility here. 
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framework, fully taking into account elastodynamic effects. The rate-and-state friction constitutive 
framework includes a single state variable </>, but is otherwise general in the sense that no special 
properties of f{4>, v) are being specified and rather generic dynamics of cj) are considered. Never¬ 
theless, we will be mostly interested in the physically relevant case in which the interface exhibits 


a positive instantaneous response to velocity changes, d.uf{4>, v)>0 (positive direct effect) (Marone 


1998 

Baumberger and Caroli 

2006), and is steady-state velocity-strengthening, dyf{4>o{v),v)>0, 

over some range of slip velocities ( 

Bar-Sinai et al. 

2014 

). In addition, we will consider two variants 


of the constitutive model, each of which incorporates a regularized response to normal stress vari- 


ations (Linker and Dieterich 

1992 

Dieterich and Linker 

19921 

Prakash and Clifton 

1992 

1993 

Prakash, 1998t Richardson and Marone, 1999 Bureau et ah, 2000) . 


While our analysis remains rather general, the main simplification we adopt is that we consider 
the limit of strong material contrast, i.e. we take one of the solids to be non-deformable. The 
motivation for this is two-fold. First, we know that the bi-material effect that emerges from the 
coupling between interfacial slip and normal stress variations becomes stronger as the material 
contrast increases (Rice et ah, 2001). We are interested here in exploring the ultimate range of 


stability and consequently we focus on strongly dissimilar materials, which will allow us to extract 
upper bounds on the stability of bi-material frictional interfaces. Second, the strong dissimilar ma¬ 
terials limit somewhat reduces the mathematical complexity involved and makes the problem more 
amenable to analytic progress, as will be shown below. We suspect that this simplification does 
not imply qualitative differences compared to the finite material contrast case, though interesting 
quantitative differences may emerge and will be explored elsewhere. The finite material contrast 
case can be studied along the same lines, though it is more technically involved. 

The structure of this paper is as follows; in Sect. the main results of the paper are listed. 
In Sect. the basic equations and the constitutive framework are introduced. In Sect. the 
linear stability spectrum for finite height H systems is derived, with a focus on standard rate- 
and-state friction. In Sect. the linear stability spectrum is analyzed in the small wave-number 
k limit, demonstrating the existence of a universal instability (independent of the details of the 
friction law) with a fastest growing mode characterized by a wave-number satisfying kH 0(1) 
and a growth rate that scales with H In Sect. the linear stability spectrum is analyzed in 
the large systems limit, kH —)> oo. We derive the stability phase diagram in terms of a relevant 
set of dimensionless parameters that quantify the competing physical effects involved. We show 
that the dominant instability mode travels at nearly the dilatational wave-speed (super-shear) 
in the opposite direction to the sliding motion. In a certain parameter range this instability is 
manifested through unstable modes at all wave-numbers, yet the frictional response is shown to be 
mathematically well-posed. Instability modes which travel at nearly the shear wave-speed in the 
sliding direction are shown to exist in a relatively small region of the parameter space. Finally, 


previous results obtained in the quasi-static regime (Rice et ah, 2001) are shown to be relevant 
within a narrow region of the parameter space. In Sect. a hnite-time regularized response to 
normal stress variations, within the framework of generalized rate-and-state friction models, is 
studied. We show that this regularized response tends to promote stability. Section offers a brief 
discussion and some concluding remarks. 


2. The main results of the paper 

The analysis to be presented below is rather extensive and somewhat mathematically involved. 
Yet, we believe that it gives rise to a number of physically signihcant and non-trivial results. In 
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order to highlight the logical structure of the analysis and its major outcomes, we list here the 
main results to be derived in detail later on: 

• The stability of a deformable body of finite height H steadily sliding on top of a rigid solid 
is studied. The linear stability spectrum is derived in the constitutive framework of velocity¬ 
strengthening rate-and-state friction models and an instantaneous response to normal stress 
variations. 

• The spectrum takes the form of a complex-variable equation implicitly relating the real wave- 
number k and the complex growth rate A of interfacial perturbations. Physically, it represents 
the balance between the perturbation of the elastodynamic shear stress at the interface and 
the perturbation of the friction stress (the latter includes the elastodynamic bi-material effect 
and constitutive effects). The spectrum can feature several distinct classes of solutions. 

• The linear stability spectrum is first analyzed in the small kH limit and the existence a uni¬ 
versal instability (independent of the details of the friction law) is analytically demonstrated. 
The instability is shown to be related to waveguide propagating modes, which are strongly 
coupled to the height H. They are characterized by a wave-number satisfying kH 0(1) 
and a growth rate that scales with H 

• As the growth rate of the waveguide-like instability vanishes in the H^oo limit, the linear 
stability spectrum is analyzed also in the large kH limit, where additional instabilities are 
sought for. Two classes of new instabilities, qualitatively different from the waveguide-like 
instability, are found: (i) A dynamic instability which is mediated by modes propagating 
at nearly the dilatational wave-speed (super-shear) in the opposite direction to the sliding 
motion and features a vanishingly small wave-number at threshold (ii) A dynamic instability 
which is mediated by modes propagating at nearly the shear wave-speed in the direction of 
sliding motion and features a finite wave-number at threshold. 

• In addition, a third type of instability — which was previously discussed in the literature — 
exists in the quasi-static regime. 

• In all cases, even when all wave-numbers become unstable in a certain parameter range, the 
frictional response is shown to be mathematically well-posed. 

• A comprehensive stability phase diagram in the large kH limit is derived, presented and 
physically rationalized. The stability phase diagram is expressed in terms of relevant set of 
dimensionless parameters that quantify the competing physical effects involved. 

• A regularized, hnite-time, response of the friction stress to normal stress variations is analyzed 
in detail and is shown to promote stability. That is, the instabilities mentioned above still 
exist, but their appearance is delayed and the range of unstable wave-numbers is reduced as 
compared to the case of an instantaneous response to normal stress variations. 

• The results may have implications for understanding the failure/rupture dynamics of a large 
class of bi-material frictional interfaces. 
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3. Basic equations and constitutive relations 


The problem we study involves an isotropic linear elastic solid, of infinite extent in the x- 
direction and height H in the y-direction, homogeneously sliding at a velocity vq in the ^-direction 
on top of a rigid (non-deformable and stationary) half space. Two-dimensional plane-strain defor¬ 
mation conditions are assumed and the geometry of the problem is sketched in Fig. 


^yy ^0 


Figure 1: A schematic 
a;-direction. 


i III 



The deformable solid is described by the isotropic Hooke’s law dij = 2^eij + {K — ‘^)6ijekk, 
where cr{x,y,t) is Cauchy’s stress tensor and the linearized strain tensor e{x,y,t) is derived from 
the displacement field u{x,y,t) according to £ij = ^ {diUj+djUi). K and y are the bulk and shear 
moduli, respectively. The bulk dynamics are determined by linear momentum balance 


djdij = pili , ( 1 ) 

where p is the mass density and superimposed dots represent partial time derivatives. A uniform 
compressive normal stress of magnitude (Jq is applied to the upper boundary of the sliding body 


ayy{x,y = H,t) = -do , (2) 

where uo > 0 and hence dyy{x, y = H,t) <0. To maintain steady sliding at a velocity uq one can 
either impose 

Uxix,y = H,t) = vo or dxy{x,y = H,t) = To (3) 

such that Vq emerges from the latter relation (for a given tq, vq is determined by the friction law, 
see below). In the limit H^oo, these two boundary conditions are equivalent as all perturbation 
modes decay far from the sliding interface. This is not the case for a hnite H, as will be discussed 
later. 

To complete the formulation of the problem, we need to specify the boundary conditions on 
the sliding interface at y = 0. First, as the lower body is assumed to be infinitely rigid, we have 


Uy{x,y = 0,t) = 0 . (4) 

Note that interfacial opening displacement, Uy{x, y = 0,t)>0, is excluded, which is fully justified in 
the context of the linearized analysis to be performed below. Equation Q is obtained in the limit 
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of large material contrast. In the opposite limit, i.e. for interfaces separating identical materials, 
the boundary condition reads ayy{x,y = 0) = —do (assuming also symmetric geometry and anti¬ 
symmetric loading). This difference in the interfacial boundary condition is responsible for all 
of the bi-material effects to be discussed below. The identical materials problem is qualitatively 
different from the bi-material problem and most of the instabilities we discuss below for the latter 
problem are expected not to exist for the former. 

Next, we should specify a friction law which takes the form of a relation between three interfacial 
quantities 

a{x,t) =-ayy{x,y = 0,t) , T{x,t) = a^y{x,y = 0,t) , u(x, t) = Ua;(x, y = 0, t) , (5) 

and possibly a small set of internal state variables (note that in the general case we have v{x, t) = 
Ux{x,y = 0~^,t) — Ux{x,y = 0~,t). Here iixix, y = 0~,t) = 0 due to the rigid substrate and only y > 0 
is of interest). As was mentioned in Sect. the class of rate-and-state friction models we consider 
involves a single internal state variable, which we denote by (i){x, t) (the physical meaning of cj) will 
be discussed later). Within this constitutive framework, the two dynamical interfacial fields, the 
friction stress r(x, t) and (f){x, t), satisfy a coupled set of ordinary differential equations in time (no 
spatial derivatives are involved) of the form 

f = F{T,a,v,4>) and cj) = G{T,a,v,(j)) . (6) 

Various explicit forms of the functions F{t, a,v, (p) and G{T,a,v,4>) will be discussed later (some 
were already alluded to in Sect. [^. 

Under homogeneous steady-state conditions, when do and vq are controlled at the upper bound¬ 
ary of the sliding body, tq and cj)o are determined from Eq. (© according to 

F{To,ao,vo,(po) = 0 and G(ro, do, uq, 0o) = 0 , (7) 

and the steady-state displacement vector field u^^\x,y,t) satisfies 

u^x\x,y,t) = vot + —y and u^y\x,y,t) = - y . (8) 

y ^ + 3 

We assume uq > 0 hereafter. Our goal in the rest of this paper would be to study the linear stability 
of this solution, for various constitutive relations. 


4. Linear stability spectrum for standard rate-and-state friction 

We consider first the standard rate-and-state friction model in which the friction stress r de¬ 


pends on the slip velocity v and a state variable cj) (Marone, 1998; Baumberger and Carol!, 2006). 


The latter is of time dimensions and represents the age (or maturity) of contact asperities (Di- 


eterich 

1978 

|Ruina 

1983 


Baumberger and Carol! 2006). It is directly related to the amount of 


real contact area at the interface (the “older” the contact, the larger the real contact area is) and 


hence to the strength of the interface (Dieterich and Kilgore 


1994 

Nakatani 

2001 

Nagata et al. 


2008 Ben-David et ah, 2010). In static situations, when u = 0, the age of a contact formed at time 


to is simply <j) = t — to. In steady sliding at a velocity vq, the lifetime of a contact of linear size 


D is (po = D/vq (Dieterich 

1978 

Teufel and Logan 

Nakatani, 2001; [Baumberger and Caroli 

2006 

). D is 


1978 Rice and Ruina 1983 Marone, 1998 
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role in this class of models. These two limiting cases are smoothly connected by choosing the 
function G{-) in Eq. ([^ such that (Ruina, 1983 Rice and Ruina, 1983 Marone, 1998; Baumberger 


and Caroli, 2006) 


(j) = G{T,a,v,(p) = g[^] , 


(9) 


with gf(l) = 0 and g{0) = 1. Furthermore, as sliding reduces the age of the contacts, we typically 
expect g\l) < 0. 

The evolution of r is assumed to take the form 


f = F(r, cr, V, (j)) = (r - v)) . 

For T—)-0, which is assumed here (in Sect, [^we will consider a hnite T), we obtain 

r = cr/(0,u) , 


( 10 ) 


( 11 ) 


which completes the dehnition of the standard rate-and-state friction model. 

Within a linear perturbation approach, the steady-state elastic helds in Eq. Q and the steady- 
state of the state variable, cpo = D/vo, are introduced with interfacial (i.e. at y = 0) perturbations 
proportional to exp [At — ikx], where A: is a real wave-number and A is a complex growth rate. The 
ultimate goal then is to find the linear stability spectrum, A{k), and in particular to understand 
under which conditions 3?(A) changes sign. 3^(A) < 0 corresponds to stability as perturbations 
decay exponentially in time, while 3?(A) > 0 corresponds to instability as perturbations grow 
exponentially. 

Before we perform this analysis, we would like to gain some physical insight into the structure 
of the problem. For that aim, we consider Eqs. © and (0 , and calculate the variation of the 
latter in the form 

^O'xy — — <^0 ^ f / ^'^yy j (1^) 

where all quantities are evaluated at the interface, y = 0, and the variation is taken relative to 
the same interfacial field (e.g. the perturbation in the slip 5ux or slip velocity 5v). When all of 
the terms are evaluated, the above expression becomes an implicit equation for the linear stability 
spectrum A(fc). R is obtained through a balance of three contributions. One contribution is 6axy, 
which is determined from the elastodynamic perturbation problem and is not directly coupled to 
the friction coefficient /. Another contribution is proportional to 5ayy, which is also determined 
from the elastodynamic perturbation problem, but which is multiplied by the friction coefficient /. 
As this is the only place in the perturbation analysis where / appears explicitly, every term that 
is proportional to / in the expressions to follow, can be physically identified with the variation of 
the normal stress ayy (and hence with the bi-material effect). Finally, the remaining contribution 
is determined by the variation of the friction coefficient 5f{(j),v), which is affected both by the 
variation of the slip velocity v and of the state variable (/>, and is multiplied by the applied normal 
stress (To. Next, we aim at calculating each of these contributions, which are being intentionally 
presented in some detail. 


4 . 1 . Perturbation of the elastodynamic fields 

We seek a solution of the linear momentum balance Eqs. Q, coupled to Hooke’s law, which 
is proportional to exp[At — ikx]. Consequently, we derive the general solution 5u{x,y,t) of the 
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problem as a superposition of shear-like and dilatational-like modes (which lead to an interfacial 
perturbation proportional to exp [At — ikx\) in the form 


6 u^{x,y,t) 

duy{x,y,t) 


-kg 

ik 


k 

-ikd 


ks 

ik 


k 

ikd 


(Ai exp[-A:sy]\ 

A2expl-kdy] 

A3 exp[/csy] 

\A4^exp[kdy] ) 


exp[At — ikx] 


(13) 


where {Aj}j=i _4 are yet undetermined amplitudes and we defined ks{A,k) = ^/A^^/c^~+~k^ and 
kd{A, k) = -y^^A2/c^~+~^. The shear and dilatational wave-speeds are Cg = \Jy/p and Cd = 
respectively. Note that kg^d are in general complex (since A is complex) and that since we adopt 
the convention that the branch-cut of the complex square root function lies along the negative real 
axis, we have 51?(A:s^rf) > 0. 

To proceed, we need to impose physically relevant boundary conditions. Recall that at y = 0 
we have 6uy{x,y = 0,t) = 0 due to the presence of an infinitely rigid substrate. We then focus 
on the case in which the velocity is controlled at y = H, i.e. Ux{x,y = H,t) = vq, cf. Eq. Q. 
Consequently, as both the normal stress ayy and tangential displacement Ux are controlled at 
y = H, the perturbation satisfies 6cryy{x,y = H,t) = 6ux{x,y = H,t) = 0. Imposing these three 
boundary conditions on the general solution in Eq. ( |13[ ) allows us to eliminate the four amplitudes 
{Aj}i=i _4 in favor of a single undetermined amplitude, and express 5axy and 6ayy at y = 0 in 
terms of 6ux- Using Hooke’s law to transform the displacement field into the interfacial stress 
components, we obtain 

dcTxy =—pkd{A,k) Gi{A,k, H) 6ux and 6ayy = i pkG2{A,k, H) 6ux , (14) 

with 

G3{A,k,H) = 


coth.{H kd) A'^/c^ 


kdkg coth{Hkd) tanh{Hkg) — 


and G2iA,k,H) = 2-^^^%^^ , (15) 


coth{Hkd) 


and recall that kg^d are both functions of A and k. This completes the elastodynamic calculation 
of 6 axy and 6 ayy at the interface. 

4.2. Perturbation of the friction law 

In the next step, we consider the perturbation of the friction law. The calculation is straight¬ 
forward, yet for completeness and full transparency we explicitly present it here. The variation of 
f{4), v) with respect to slip velocity perturbations takes the form 


< 5 / = {dxf + 1 


(16) 


where all of the derivatives are evaluated at the steady-state values corresponding to vq. As the 
perturbation in (j) takes the form (5(/)~exp[At — ikx], we can use Eq. Q to obtain 


6 v 


9'(1) 


|9'(1)I 


Vo (A - y'(l)^ 


vo{A + \g'{l)\'^ 


(17) 
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where in the latter we used g''(l) = —1^'(1)| because g''(l)<0. Finally, since dv(j)o = —D/vq, we can 
relate d^f to dyf according to 

dvf = dj - . (18) 

% 


Substituting Eqs. 0 -(|18[) into Eq. ( |16[ ) and using 5v = ASux, we obtain the final expression for 
the perturbation of v) 


df = 


A 


A + 


b'(i)ko 


dyf A dyf ) 5ux 


(19) 


D 


4-3. The linear stability spectrum and dimensionless parameters 

Now that we have calculated the perturbations of 5axy, dayy and 6f in terms of 6ux, we are 


ready to derive the linear stability spectrum. Eor that aim, we substitute Eqs. (14) and (19) into 


Eq. (12) and eliminate 5ux to obtain 


SiA,k,H) = f,kd{A,k)Gi{A,k,H)-iiakfG2 (A, k,H)+ = 0 • 

A H \ / 

( 20 ) 

This is an implicit expression for A{k), which is in general a multi-valued function with various 


branches. Analyzing Eq. (20) is the major goal of this paper. 


As a hrst step, we briefly discuss the symmetry properties of Eq. (20). The latter satisfies 


S{A, —k, H) = S{A,k, H), where a bar denotes complex conjugation. This implies that for each 
mode with A:<0, there exists another mode with fc>0, which has the same growth rate 51?(A) and 
an opposite sign frequency —^(A). Consequently, we have A{—k) = A{k). These two modes have 
the same phase velocity A{A)/k (because both the numerator and denominator change sign when 
k^—k) and they propagate in the same direction. These symmetry properties allow us to assume 
hereafter A: > 0 without loss of generality. 

We would now like to discuss the set of dimensionless parameters that control the linear stability 
problem. In addition to / itself, which is obviously a relevant dimensionless quantity, we dehne 
the following three quantities 


A = ^ 

dyf 


/? = - 
Cd 


I' 


T 


Cs CTO dyf 


( 21 ) 


which involve material/interfacial parameters, interfacial constitutive functions that may depend on 
the sliding velocity (i.e. dyf and dyf), and the applied normal stress do. The first quantity. A, is the 
ratio between the steady-state and the instantaneous variation of / with v (“instantaneous” here 
means that the slip velocity variation takes place on a time scale over which the state variable (p does 
not change appreciably). As such, A is a property of the frictional interface. As mentioned above, 
frictional interfaces generally exhibit a positive instantaneous response to slip velocity changes (a 
positive “direct effect”), dyf>0, which is assumed here. Consequently, the sign of A is determined 
by the sign of dyf, that is, by whether friction is steady-state velocity-weakening or velocity¬ 
strengthening in a certain range of sliding velocities. While steady-state velocity-weakening is 


prevalent at small sliding velocities and is very important for rapid slip nucleation (Rice and 
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Ruina, 1983 Dieterich, 1992; Ben-Zion 


state velocity-strengthening (Noda and Shimamoto 


2008), some frictional interfaces are intrinsically steady- 

Moreover, it 


2009 

Ikari et al. 

2009 

2013 


has been recently argued (Bar-Sinai et al. 2014) that dry frictional interfaces generically exhibit 


a crossover from steady-state velocity-weakening, dyf < 0 , to velocity-strengthening, dyf > 0 , with 
increasing v beyond a local minimum. This claim has been supported by a rather extensive set of 


experimental data for a range of materials (Bar-Sinai et al., 2014) 


From this perspective, it would be instructive to write A using Eq. (18) as 

dy 


A = 1 - 


Dd4>f 

vldj 


= 1 - 


logfli 


/ 


d. f 

logi; ^ 


( 22 ) 


where generically d^f>0 (i.e. frictional interfaces are stronger the older - the more mature - the 
contact is) and recall that all derivatives are evaluated at steady-state corresponding to a sliding 
velocity vq (i.e. v = vq and ipo = D/vq). Hence, a crossover from A < 0 at relatively small uq’s 
to A > 0 at higher vq’s implies a crossover from to with increasing 

Vq. While our analysis can be in principle applied to any A, the most interesting physical regime 
corresponds to A>0, i.e. to sliding on a steady-state velocity-strengthening friction branch, where 
friction appears to be stabilizing and the destabilization emerges from elastodynamic bi-material 
and finite size effects (see below). For steady-state velocity-weakening friction, which has been 


studied quite extensively in the past (Dieterich, 1978 Rice and Ruina, 1983 Marone 1998), we 


have A<0 and some unstable modes always exist. This instability is of a different physical origin 


compared to the instabilities to be discussed below. Physical considerations (Bar-Sinai et al., 2014) 


indicate that increasing the steady-state sliding velocity vq along the velocity-strengthening branch 
is accompanied by either a decrease in or an increase in or both. In the limiting case, 

^log,#,/^ have A—>- 1 . Consequently, on the steady-state velocity-strengthening branch 

we have 0 < A < 1. 

The second quantity, /?, is the ratio between the shear wave-speed c* and the dilatational wave- 
speed Cd and hence is a purely linear elastic (bulk) quantity. (3 is only a function of Poisson’s ratio 
V (in terms of the shear and bulk moduli, we have v= which for plane-strain conditions 

takes the form 


/? = 


I l-2u 

2(1-i^) 


(23) 


For ordinary materials we have 0 < ^ (recall that thermodynamics imposes a broader constraint, 

— 1 < z/ < ^, but negative Poisson’s ratio materials are excluded from the discussion here), which 
translates into 0</3<^. 

The third quantity, 7 , is the ratio between the elastodynamic quantity fi/cg — proportional to 
the so-called radiation damping factor for sliding (Rice, 1993; Rice et ahj 2001 ;[&upi and Bizzarri 


2013) — and the instantaneous response of the frictional stress to variations in the sliding velocity, 


dyT = aQdyf. The latter is a product of the externally applied normal stress (Tq and dyf. As such, 
7 quantihes the relative importance of elastodynamics, the applied normal stress and the direct 
effect (an intrinsic interfacial property). It has been shown that many frictional interfaces are 
characterized by an instantaneous linear dependence of / on log v over some range of slip velocities 


9, f 

logtz •' 


IS a 


due to thermally-activated rheology. In this case, is a positive constant and dyf- 

decreasing function of vq. Therefore, as vq increases (for a fixed uo), elastodynamics becomes more 
important and 7 increases. Finally, note that as cto> 0 and dyf>0, we have 7 > 0 . 
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In terms of the four independent dimensionless parameters A, (3, 7 and /, the linear stability 
spectrum in Eq. (20) can be rewritten as 


S{A, k,H) = kd{A, k) Gi(A, /c, H) — i iik f G 2 (A, k,H) + iJi 


A A(A) + A 
Cs 7 (A(A) + 1) 


= 0 


(24) 


where A (A) = The shear modulus /r clearly drops as a common factor and then the only 

dimensional quantities left are k and A/cg, which obviously form a dimensionless combination. 

The main goal of most of the remaining parts of this paper is to find the solutions (i.e. roots) 
of Eq. (24), especially those with 3f?(A) > 0. Equation (24) is a complex equation (both in the 


complex-variable sense and in the literal sense) which includes branch-cuts in the complex-plane 
and is expected to have several solutions A{k). In looking for these solutions, one can follow various 
strategies. One strategy would be to numerically search for these solutions. We do not adopt this 
strategy in this paper. Rather, we will show that solutions of Eq. (24) are physically related to 


various elastodynamic solutions and that this insight can be used to obtain a variety of analytic 
results (which are then supported numerically). Below we treat separately the small and large kH 


limits of Eq. (24). 


5. Analysis of the spectrum in the small kH limit 


Our first goal is to analyze the linear stability spectrum in Eq. (24) in the small kH limit, where 
kH^O{l) or smaller. In this range of wave-numbers k, perturbations are strongly coupled to the 
finite boundary at y = H. The most notable physical implication of this is that the elastodynamic 
solutions in the bulk are non-decaying in the y-direction. Consequently, we will look for solutions 
that are related to waveguide-like solutions featuring a propagative wave nature in the x-direction 
and a standing wave nature in the y-direction. 

In general, a 2D wave equation would give rise to a dispersion relation of the form A^ = 
—c^(/c^ -|- ky) (c is some wave-speed). In the presence of a finite boundary at y = H, kx ^ k is 
continuous and ky is quantized according to ky = rmr/H, with m being a set of integers/half- 
integers which are determined from the boundary conditions. This quantization implies that in the 
long wavelength limit /c—>-0, the dispersion relation results in a cutoff frequency A = Aicrmr/H. 
This property will be shown below to have signihcant implications on the stability problem. 

Based on the idea that waveguide-like solutions might be important for the stability problem, 
we seek now solutions to Eq. (24) in the limit A; —)• 0. Eor that aim, we look first for solutions 
corresponding to 6axy{A, k) = 0 , which defines the relevant waveguide dispersion relation A^g(fe). 


Eollowing Eqs. (14)-(15), solutions of interest correspond to coth.{Hkd) = 0, i.e. 


Hkdi^Axug, k') — Hi 


A2 

' -^^wg 


+ F = ± 


i{2n + l)vr 


Ao = Auigik^O) = ± 


i{2n + l)vrcrf 


2H 


where Aq is the waveguide cutoff frequency. 


(25) 


Our strategy will be to expand ^(A, k, H) of Eq. (24), as a function of the two variables k and 


A, to leading order around (A = Ao,A; = 0). That is, we expand S{A,k,H) to linear order in k 
and in 5A, A~Ao -|-(^A-|- 0(<5A^), and then set S{A,k,H) = 0 to obtain 5A{k). This procedure is 
expected to yield a solution that satisfies 6A{k)^k. In principle, if indeed a solution that satisfies 
3^[5A(A:)] ~ A: exists (i.e. k = 0 is a regular point where a linear expansion exists), then irrespective 
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of the sign of the proportionality coefficient we have 3^[JA(/c)] >0 for either k>0 oi k <0, implying 
an instability. In fact, a similar conclusion can be reached even if the discussion is restricted to 
A:>0. In this case, if A(/c) ~ Ao+(fA(A;) — with a pure imaginary Aq and a complex 6A{k)^k — is 
a solution of S{A{k), k, H)=0 then the symmetry property ^(A, —k, H) = S{A, k, H) implies that 
5(Ao + I?[JA(A:)] — f3'[(5A(A:)], —k, H) = 0. Therefore, 

A{k) ~ Ao + 5R[<5A(-A:)] - i^[5A{-k)] = Aq - 3?[(5A(A:)] + i^[5A{k)] (26) 


is also a solution of S{A{k), k, H) =0. 

We thus conclude, based on the last argument, that if a solution with a complex 5A{k) ~ k 
near Aq (purely imaginary) exists, then actually there exist two solutions with the same A'[5A(A:)] 
and 51?[5A(/c)] of opposite signs. Both solutions propagate with the same group velocity d^[6A]/dk, 
but one is stable (3?[(5A(A:)] < 0) and the other is unstable (iR[(5A(/i;)] >0). This leads to the quite 
remarkable conclusion that steady-state sliding along (strongly) bi-material frictional interfaces in 
this broad class of constitutive models is universally unstable. 

To fully establish this important result, we derive it by an explicit calculation. To that aim, as 


explained above, we need to expand S{A, k, H) of Eq. (24) to leading order around (A = Ao, A: = 0). 
This is done in detail in Appendix [Appendix A[ but the essence is given here. The contribution 
related to 5ayy in Eq. ( pd] ) is proportional to both G 2 {A,k, H) and k\ the former diverges in the 
limit {A^k) —)• (Ao,0), while the latter clearly vanishes. In total, the term proportional to dayy 
approaches the finite limit 


ik f G2(A, k,H) 


fcsk 


fi'^H tan(fl| Ao|/cs)(5A 


(27) 


where the ± corresponds to Aq = ±i|Ao|. This leads to the non-trivial situation in which the 
leading order contribution involves the ratio of 5A and k. 

Consequently, the contributions related to duxy and 5f can be taken, albeit with some care, to 


zeroth order, yielding (see Appendix Appendix A for details) 


l^ol 

Cs tan(77|Ao|/cs) 

^0 Aq- l-A |Ao| / — IAo| (1 — A) ± f(A-|-|AoP) 


kd{A,k) Gi{A,k,H) 
ao6f ~ 


Cs 7(Ao + l) Cs'y 


l + |Ao|2 


(28) 

(29) 


where Aq = \g^i)\vQ ■ Collecting all three contributions we end up with 


S{A,k,H) ~ 


|Ao| 


± 


fcgk 


Aq Aq -|- a 

Cstan{H\Ao\/cs) ^ P'^Htan{H\Ao\/cs)SA ' 7(Ao- h l) 




= 0 . 


(30) 


Equation (30) clearly establishes a linear relation between 6A and k. Extracting A(A:) = Ao -|- 


6A{k) by solving Eq. (30), we obtain 


iR(A) ~ T 


fcsk 

\ |Ao| 

Aq Ao (1 - A) 


P‘^Htan{H\Ao\/cs) 

Cs tan{H\Ao\/cs) 

Cs 7 l-b|Aop 



|Ao| 


_ |Ao| |Ao|(l - A) 

Cstan(i7|Ao|/c^) l-F|Aop 
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[AqI a -I- |Ao|" 

Cs7 1 -b |AoP 


+ 0{k^ 


(31) 




























where ^[JA] can be easily obtained as well and is independent of the sign of A[Ao]. Equation (31) 
has precisely the predicted structure, i.e. there are two solutions for A in the small A: > 0 limit, 
whose real parts have opposite signs. Consequently, a solution with 3f?[A] >0 always exists, i.e. the 
system is universally unstable. 

While Eq. (31) can be somewhat simplihed, it is retained in this form so that the physical 
origin of the various terms will remain transparent. Note that in the particular case of /3 = 1/2 
(which will be used in some of the numerical calculations below), a significant simplification is 
obtained, leading to an unstable branch 3f?(A)~4//cCs/[(2n + l)7r] + 0{k'^) (some care should be 
taken when obtaining this result as a naive substitution of = lj2 in Eq. (31) results in some 
apparently divergent contributions). The analytic result for the small k behavior of the growth 
rate il?[A] presented in Eq. ([^, which is one of the main results of this paper, is verihed in Eig. [| 
for the few first n’s by a direct numerical solution of the linear stability spectrum in Eq. (24). The 
numerical solution of the spectrum shows that the most unstable mode satisfies kH^0(1), where 
iR[A] attains its first maximum (corresponding to the n = 0 solution). This can be analytically 
obtained by calculating the 0{k‘^) correction to Eq. (31), though the calculation is lengthy. 



Figure 2: (left) The growth rate 5R[A] (in units of Cs/H) vs. Hk, for small Hk, for various values of 7 / (solid lines) 
and quantization numbers n (dashed and dotted lines). 7 / was varied by varying 7 and we used H = 0.5, where H 
is H measured in units of CsD/\g'{l)\vo, which is the product of the basic time scale of the system, D/\g' {1)\vq, and 
the velocity Cs. The black broken lines show the asymptotic behavior at fc—>0, as predicted by Eq. (311. (right) The 
same as the left panel, but for varying H (and n = 0). Note the saturation for large H, as predicted analytically. 
Unless specified otherwise, all calculations, here and in what follows, were done with the generic parameter values 
f = 0.3, 7 = 3, /3 = 0.5 and A = 0.7. 


We thus conclude that the finite size H of the sliding system has significant implications for 
its stability, in particular it implies the existence of an instability with a wavelength determined 
by H. This instability should be relevant to a broad range of systems, for example an elastic 
brake pad sliding over a much stiffer substrate, for which recent numerical results demonstrated 
dominant instability modes directly related to the intrinsic vibrational modes of the pad (]Behrendt 


et ah, 2011 Meziane et ah, 2007). The universal existence of this finite H instability does not 


immediately mean that it will be indeed observed since other instabilities, which do not necessarily 
satisfy kH^0(1), might exist and feature a larger growth rate (when several instabilities coexist, 
the one with the largest growth rate will be the dominant one). 

To further clarify this point, we consider large, but finite, H in Eq. (31). 


In this limit, by 


counting powers of H and substituting k ^ H ^ for the fastest growing mode, we obtain for the 
latter ^[A] oc f Cg/H. This scaling is verified numerically in Fig. (right panel). This result shows 
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that this instability depends on the presence of friction, but not very much on the details of the 
friction law (e.g. the length scale D does not play a dominant role), and that the growth rate of 
the instability decreases with increasing H. This raises the question of whether in the large H 
limit there exist other instabilities with a larger growth rates, which will be addressed in the rest 
of the paper. 

It is important to note that while we highlight the role of the hnite size H in relation to 
the universal instability encapsulated in the growth rate in Eq. (31), we should stress that the 


bi-material effect remains an essential physical ingredient driving this instability. This is evident 
from the observation that 5?[A] is proportional to / in Eq. (31), where — as explained in Sect.— 


the latter is a clear signature of the variation of the normal stress with slip, which is associated with 
the bi-material effect. We thus expect that this universal instability does not exist for frictional 
interfaces separating identical materials. The effect of finite material contrast should be assessed 
in the future. 

Before we discuss the large H limit, kH —oo, we would like to note another interesting 
implication of the finite system size H. Eor infinite systems, H ^ oo, there is an equivalence 
between velocity-controlled and stress-controlled external boundary conditions (cf. Eq. ([^) because 
perturbations decay exponentially away from the interface in the y-direction. This equivalence 
breaks down for a finite H. The analysis above focussed on velocity-controlled boundary conditions. 
In Appendix Appendix B, we consider also stress-controlled boundary conditions and explicitly 


demonstrate the inequivalence of the two types of boundary conditions for finite size systems. The 
differences, though, are quantitative in nature and the generic instabilities discussed above remain 
qualitatively unchanged. 


Einally, we would like to note that there can be solutions to Eq. (24) in the small kH limit 


other than Eq. (31). We are not looking for them here because the result in Eq. (31) already shows 


that the system is always unstable. In addition, we expect the decay of the growth rate 3^[A] with 


i/ to be a generic property of unstable solutions of Eq. (24) in the small kH limit. Consequently. 


we focus next on the large kH limit, looking for qualitatively different unstable solutions. 


6. Analysis of the spectrum in the kH oc limit 


After analyzing the linear stability spectrum of Eq. (24) in the small kH limit, our goal now is 
to provide a thorough analysis of the opposite limit, kH —)-oo. The length H enters the problem 


through the elasticity relations in Eq. (14) and more precisely through the functions Gi ^2 in Eq- (15). 
Taking the kH^oo limit in Eq. (15), which amounts to taking the arguments of coth(-) and tanh(-) 
to be arbitrarily large, we obtain 


Gi(A, k, H^oo) 
G2(A, k, H^oo) 


A2 


gi{A,k) ^ksiA, k) kd{A, k) - k^] ' 
92 {A,k) = 2- gi{A,k) , 


(32) 


which should be used in Eq. (24). The friction part is of course independent of H. 


To further simplify the analysis of the spectrum in this limit, we define an auxiliary (and 
dimensionless) complex variable z that relates the spatial and temporal properties of perturbations 
according to 
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Defining the dimensionless wave-number as g = , and recalling that we already defined above 

the dimensionless (complex) growth rate as A = 
these definitions, Eq. (24) can be rewritten as 


9'(i)ho ’ 


Eq. (33) can be cast as X = —iqz. With 


'S(^, g) = 7 (1 - iQz) - /3'^z‘^ gi{z, /3) - ifg 2 {z, /3)j - iz (A - iqz) = 0 , 


where 


(34) 


9i{z,I3) = 


1 — yjl — z‘^\/l — I3‘^z‘^ 


and g2{z,l3) = 2-gi{z,P) . 


(35) 


As before, Eq. (34) is an implicit expression for the explicit spectrum z{q)^ which depends on 
the four dimensionless parameters A, /3, 7 and /. Note also that due to algebraic manipulations. 


Eq. (34) no more follows the structure of Eq. (12) (which is preserved in Eq. (20) and (24)), rather 


the terms are mixed to some extent. Later on, when discussing some of the physics behind our 


results, we will reinterpret them in terms of Eq. (12). Due to the appearance of the complex 


square root function in the above expressions, Eq. (34) is understood as having a branch-cut on 


the real axis along | 2 :| > 1 (there is also a branch-cut on the real axis along \z\>l/ fi associated with 

— P'^z'^. Combinations of \/l — and y^l — P^z"^, as in Eq. ( [Sd] ), may have more complicated 
branch-cut structures). The existence of these branch-cuts has implications that will be discussed 
later. Finally, in analyzing Eq. (34) we assume that the dimensionless wave-number q — and hence 
the dimensional wave-number k — spans the whole interval 0 < g < 00 . The small wave-numbers 
limit, g—tO, is understood to imply Dk<^\g'(1 )\vq/cs while maintaining kH^l. This can always 
be guaranteed by having a sufficiently large H. 

In the next parts of this section we present an extensive analysis of the linear stability spectrum 
in Eq. (34). As in Sect. we will establish relations between the unstable solutions of Eq. (34) 
and various elastodynamic solutions and use this insight to derive analytic results that will shed 


light on the underlying physics. In Sect. 6.1 we show that there exist unstable solutions related to 


dilatational waves propagating in the direction opposite to the sliding motion. In Sect. 6.2 we show 
that there exist another class of unstable solutions which are related to shear waves propagating 


in the direction sliding motion. In Sect. 6.3 we briefly review a qualitatively different class of 


solutions, which are not elastodynamic in nature, but rather quasi-static (Rice et al., 2001). In 


Sect. 6.4 we present a comprehensive stability phase diagram, which puts together all three classes 


of solutions of the linear stability spectrum in Eq. (34). While we do not provide a mathematical 


proof that other classes of solutions do not exist, we suspect that our analysis is exhaustive. 

6.1. Dilatational wave dominated instability 

In the spirit of Sect. we will look for solutions of Eq. (34) that are related to propagating 
wave solutions. In particular, we note that the linear stability spectrum of Eq. (34) significantly 


simplifies when 2 ; = 1//3, for which y^l — vanishes. Physically, the latter corresponds to 
dcTxy = 0, where kd{A,k) = 0 and gi{K,k) is finite (cf. Eq. ([IH ) with gi replacing Gi), i.e. to 
frictionless boundary conditions. Substituting z = \/l3 in Eq. ( |34[ ) and taking the limit g—tO, we 
immediately observe that it is a solution if 7 /(l// 3 ^ — 2 ) = A//3. As will be shown soon, the latter 
is an exact stability condition for the emergence of unstable solutions located near z = \lfi in the 
complex z-plane. Recall that a real 2 ; is equivalent to 3^(A) = 0, which is precisely where solutions 
change from growing (iR(A)>0) to decaying (3f?(A)<0) in time. 
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z = \l(5 corresponds to the dispersion relation for dilatational waves, A = —ic^k (i.e. a frictionless 
boundary conditions, 5axy = 0), which means that instability modes located near z = l//3 in the 
complex-plane travel at nearly the dilatational wave-speed in the direction opposite to the sliding 
direction. The direction of propagation is a result of the minus sign in the last expression. It is 
important to stress in this context that while 2 ; —)■ —1//3 also corresponds to a dilatational wave 


{Scrxy = 0), it is not a solution of Eq. (34). That is, friction in the presence of homogeneous sliding 
breaks the directional symmetry of dilatational waves. The fact that A vanishes in the limit fc—)-0, 
marks a crucial difference between the analysis to be performed here and the one in Sect.[^ where 
the finite system size H implied a hnite cutoff frequency A—t-Aq in the limit /c—)-0. 

Following this simplified analysis, which indicates that some unstable solutions might be located 
near z = l//3 in the complex z-plane, we aim at obtaining analytic results for the spectrum by a 
systematic expansion around this point. That is, we are interested in obtaining a systematic 
expansion of the form z = 1/(3 + 6z, where 5z is a small complex number (i.e. \5z\ <C 1) whose 
imaginary part determines the stability of sliding (A(<5z) > 0 implies instability and '3s{5z) < 0 
implies stability). This should be done carefully, though, since z = l/l3 is a branch-point, where a 
Laurent expansion does not exist. 

To address this issue, let us briefly discuss one of the physical implications of being close to 
z = l/(3. First, note that the real part of kd in the elastodynamic solution in Eqs. (13) controls 


the decay length in the y-direction (kg plays a similar role, but is not discussed here. Note also 
that ^ 3^4 —)• 0 in the H ^ 00 limit considered here, which ensures the proper decay of solutions 
sufficiently away from the interface.). Then, expressing kd in terms of z, kd = ky^l — (3^z‘^, we 
observe that kd vanishes as z^l/f3, i.e. there is no decay in the y-direction in this case, and the 
proximity to z = l/(3 actually controls the smallness of kd (for a given wave-number k). Therefore, 
we define a complex number Kd^ \Jl — (3‘^z‘^ such that kd = k Kd, where <C 1. 

With this definition of smallness, we go back to our original motivation to derive a systematic 
expansion around z = l/(3 and express z in terms of Kd as 


z = 




/3 


1//3 - ^ + 0 ( 4 ) 


(36) 


The latter expression has the desired form z = l/(3 + Sz and our next goal is to estimate Kd itself 
from the linear stability spectrum in Eq. (34). To do this, we need to rewrite the spectrum in 


terms of the new independent variable Kd- In the proximity of z = l/(3, the functions gi{z) and 
g 2 {z) in Eq. (35) can be written in terms of Kd as follows 




1 /( 3 ^ 


1 =F ^Kd^/l/W^ 


and 




(37) 


where the minus sign corresponds to the stable branch (A( 2 ;) < 0) and the plus sign to the unstable 
branch (A( 2 ;) >0). This emerges from the limit —Vl — —)• — 1 as 1/(3, where the 

different signs correspond to taking the limit from the two sides of the branch-cut (^(^z) —)• 0^). 
The main advantage of Eqs. (|37| ) is that gi^{Kd) and g 2 ^{>^d) are analytic such that a Laurent 
expansion around Kd = 0 exists 


^We note in passing that we could do the whole analysis with Sz instead of Kd, invoking the leading term ^y/Sz 
in a fractional power series. The two routes are equivalent when identifying Kd — ■\/—2g~Sz, which is precisely what 
Eq. (361 states. 
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Using Eqs. (37), we can rewrite the linear stability spectrum in Eq. ( |34[ ) in terms of as 

siKd,q) 7(1 - iQ/13) [ndgi^ii^d) - if h^ind)] - i/PiA - iq/P) = 0 , ( 38 ) 

where we set z = ljl3. We then linearize the following K^-dependent quantities 

kdgi^iKd)czKd/l3‘^ + 0 {kI), g2^{Kd)zz2-l/l3'^ (l±iKdy/l/W^'^+0 {kI) , (39) 

substitute them into Eq. (38) and solve the resulting linear equation for Kd, obtaining 

i /7 (1 - iq/P) (2 - l//3^) + i (A - iq/P)/P 
7 (1 - V/3) (1 T/vW^)//?' 

Substituting the latter in Eq. ( p^ , we can calculate the dimensionless growth rate iR(A) = q^{z) 
in the form 


q"^ {1 - A) 


' [7//3 (l//j^ - 2) - l]gV/3^ + 7//3(1//3^ - 2) - 
72 (9V/32 + l)2(l^/^l/;32-l)2 


(41) 


where, as before, the stable solution corresponds to the minus sign and unstable one to the plus 
sign. This analytic prediction is one of the major results of this paper. It is important to stress 
that unlike the growth rate in Eq. (31), which was obtained by a small wave-numbers expansion, 
the growth rate in Eq. (41) was obtained by an expansion in the complex plane near z = 1//3. 
Consequently, it is valid — as will be explicitly demonstrated below — for any wave-number q. 

A lot of analytic insight can be gained from Eq. (41). First, note that the growth rate iR(A) in 
Eq. (41) is continuous, but not differentiable at the transition from the stable to unstable branches 
as a function of q (i.e. it has a kink due to the existence of a branch-cut in the equation for the 
spectrum). Then, we see that unstable modes appear in the long wavelength regime, 0<q<q^'^'>, 
where the critical wave-number q^‘^'> is simply obtained from the condition I?(A) =0 


q(d) ~ 


hf/3{l//3^-2)-A 

1 - 7 // 3 ( 1 // 32 - 2 ) 


(42) 


The instability threshold is obtained by taking the limit q^‘^'> —^0 (i.e. the instability does not occur 
at a finite wavelength) 

^ /3(1//32-2) ’ 

which is identical to the one derived at the beginning of this section. When the left-hand-side is 
smaller than the right-hand-side, there exist no unstable modes, i.e. the regime (XqKq^f) shrinks 
to zero and iR(A) is always negative. In the opposite case, when the left-hand-side is larger than 
the right-hand-side, a finite range of unstable modes emerges. A simple calculation shows that the 
threshold condition actually emerges from the numerator of in Eq. (40), and in particular 

from its g-independent part (since the threshold condition corresponds to the limit of vanishing 
wave-number q). 

Let us discuss the physics embodied in Eq. (43). For that aim, recall the definitions of 7 and 


A in Eq. (21) and substitute them in Eq. ( |43[ ) to obtain 

/^/3(1//32-2) =c7od,/ . 

Cs 
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A first observation is that this instability threshold is independent of dyf. Put differently, as far as 
the threshold is concerned, the distinction between dyf and dyf is irrelevant as if the friction law 
is only rate-dependent, f{v). This can be understood as follows; the right-hand-side of Eq. (44) 


corresponds to the cto 5/ term (variation of the friction law) in Eq. (12). Near threshold we have 
A~—i/c—)-0, which can be substituted in the expression for 6f in Eq. (19). We observe that the 


term proportional to dyf scales as k'^, while the one proportional to dyf scales as k and hence the 
latter dominates the former. Consequently, we have Sf^dyf, independently of dyf. 


Another aspect of Eq. (44) which is worth noting concerns the left-hand-side, which is pro¬ 
portional to / and hence corresponds to the fdayy term in Eq. (12). Indeed, 5ayy in Eq. (14) 


scales as k, while 5axy is higher order in k and hence negligible. Moreover, we observe that 5a. 


yy 


is proportional to the so-called radiation damping factor for sliding (Rice, 1993; Rice et ah, 2001 


Crupi and Bizzarri, 2013| ), fi/cg, which essentially follows from dimensional considerations in the 
elastodynamic regime. To conclude, the present discussion shows that Eq. (43) is actually of the 
form aodyf ^ f 5ayy, i.e. the onset of instability is controlled by a balance between the stabilizing 
steady-state velocity-strengthening friction, dyf, and the destabilizing elastodynamic bi-material 
effect, 5ayy. 

Next, we consider the analytic prediction for 3f?(A) in the limit k ^ oo (or in dimensionless 
units, 3?(A) in the limit q —)• oo). By counting powers in Eq. (41) we observe that in the limit 
(?—)-oo, 3?(A) approaches a constant whose sign is determined by the sign of jf/3 (l//3^ — 2) — 1. 
If ^fl3 (l//3^ — 2) < 1, then the constant is negative. This does not immediately imply stabil¬ 
ity because, following the discussion above, a finite range of unstable modes emerges if A < 
7//3(1//32-2)<1. If, on the other hand, we have 


7//3 (1//3' - 2) > 1 , 


(45) 


then the constant is positive, which implies that all wave-numbers are unstable. Indeed, Eq. (42) 
shows that the critical wave-number diverges in the limit ■yflS (l//3^ — 2) —)-l. 

We have thus seen that when elastodynamic effects become sufficiently strong, i.e. when the 


combination 7/ becomes sufficiently large, all wave-numbers are unstable. This observation raises 


the issue of ill-posedness, which has been quite extensively discussed in the literature recently (Re¬ 


nardyl 1992[ 

Adams, 1995[ Martins and Sim6es[ 1995; Martins et al., 

19951 Sinioes and Martins, 

1998[ Ranjith and Rice 

2001 

). Ill-posedness is a stronger condition than instability for all wave- 


numbers, i.e. a problem can feature unstable modes at all wave-numbers but still be mathematically 
well-posed, and is defined as follows; consider the perturbation of any relevant interfacial field in 
the linear stability problem, e.g. the slip velocity field 5v, and express it as an integral over all 
wave-numbers 

POO 

5v{x,t) ^ / a{k) exp[—ikx] exp[A{k)t]dk , (46) 


where a{k) is the amplitude of the fcth mode. If this integral fails to converge, the problem is 
regarded as mathematically ill-posed. 


An important example in this context (Renardyl 

19921 |Adams[ |1995t |Martins and Simoes 

1995 

Martins et al. 

1995 Simoes and Martins 

1998 

Ranjith and Rice, 12001) is sliding along a 


bi-material interface described by Coulomb friction, T = af (where / is a constant). In this case, 
3?(A(/c)) ~ \k\ (with a positive prefactor) and the integral in Eq. (46) fails to converge for any 
X 7^ 0 at any finite time, unless a{k) decays exponentially or stronger with |fe|. The problem can 
be made well-posed if in response to normal stress variations, r = af is approached over a finite 
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time scale (Ranjith and Rice, 2001). In our problem, within the standard rate-and-state friction 
framework, we saw above that there exists a range of parameters in which all wave-numbers are 
unstable. Yet, in this case ^{A{k)) approaches a constant as k^oo, in which case the integral in 
Eq. (46) converges. Therefore, we conclude that the response of bi-material interfaces described 
by standard rate-and-state friction laws is mathematically well-posed. 

We are now in a position to quantitatively compare the analytic predictions derived from Eq. 
(41) to a direct numerical solution of the linear stability spectrum in Eq. (34). The results are 
shown in Eig. On the left panel, 51?(A) = q^{z) is shown as a function of q for various y/’s and 
fixed representative values of A, / and (3. 



Figure 3: (left) The dimensionless growth rate 5R(A) = qA{z) vs. the dimensionless wave-number q for various y/’s 
( 7 / was changed by changing 7 ). The solid lines show the numerical solution of the linear stability spectrum in 


Eq. (341 for both the stable (K(A) < 0) and unstable (5R(A) > 0, when it exists) branches. The dotted lines correspond 


to the analytical prediction of Eq. (411. A very good quantitative agreement between the analytic prediction and 


the direct numerical solution is demonstrated (see text for more details). The discontinuities (gaps) observed in the 
full numerical solution are discussed in Appendix Appendix C (right) 3 ff( 2 ) vs. the wave-number q. The parameters 
used are f = 0.3, (1 = 0.5 and A = 0.7. 

All in all, Fig. demonstrates a good quantitative agreement between the analytic prediction 
and the full numerical solution over a significant range of parameters and wave-numbers. In 


particular, Eq. (43) predicts (for the chosen (3) that the onset of instability takes place at 7/= 0.7, 
which is precisely what is observed. Furthermore, the onset of instability appears at fe —)• 0, as 
predicted. The critical wave-number q^<^'>, predicted in Eq. (42), is quantitatively verified for 


several sets of parameters above threshold. Finally, the shape of the unstable spectrum, including 
the constant asymptote as g—)-oo, is quantitatively verified. 

The only interesting deviation of the analytic prediction of Eq. (41) from the full numerical 


solution in the left panel of Fig. is that the latter exhibits a discontinuity (a gap) at the transition 
from the unstable to the stable part of the solution, while the former is continuous but rather 
exhibits a discontinuous derivative. This results in a shift of the stable part of the spectrum when 
an unstable range of wave-numbers exists. The origin of the gap in the spectrum is explained in 
Appendix Appendix C On the right panel of Fig. 3?(2;) is shown as a function of q for both the 
numerical solution of Eq. (34) and the real part of Eq. (|36[) (together with Eq. (40)). The figure 


demonstrates, again, a good quantitative agreement between the analytic prediction and the exact 
numerical solution. Furthermore, the two panels of Fig. [^show that indeed 9‘(2:) = 0?(A)/g^)l?(2;)~ 
1//3, as expected for solutions located near z = l//3. In particular, note that solutions remain close 
to z = \/(3 in the complex-plane for every wave-number in this class of solutions. 
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With this we complete the discussion of the dilatational wave dominated instability, which cor¬ 
responds to unstable modes of predominantly dilatational wave nature propagating in the direction 
opposite to the sliding direction (corresponding to solutions near 2 : = l//3). In the next subsection 


we discuss a distinct class of unstable solutions of the linear stability spectrum in Eq. (34). 


6.2. Shear wave dominated instability 

Inspired by the discussion in the previous subsection, we look here for another class of unstable 
solutions. This time we focus on the zeros of \/l — z'^, in particular on solutions located near z = — l 
in the complex-plane. As will be shown below, these instability modes are of shear wave-like nature, 
propagating with a phase velocity close to Cg in the sliding direction. 


To see how this rigorously emerges, we set z = — l in Eq. (34) (which corresponds to 51?(A) = 0 


i.e. to the threshold of instability) and separate the real and imaginary parts to obtain 


= 


7v/r^ 

1-7/ 


and 


A = 7/ - 




(47) 


where g7) is the critical (dimensionless) wave-number at threshold and the second relation is the 
onset of instability condition (an instability occurs when the left-hand-side is smaller than the 
right-hand-side). This is an exact result. Unlike the dilatational wave dominated instability, 
which featured a vanishing critical wave-number at threshold, —>- 0 , the shear wave dominated 
instability takes place at a finite wave-number (above threshold, a finite range of unstable q’s 
emerges around cf. Fig. |^. 2 : = —1 corresponds to the dispersion relation for shear waves, 
A = icsk, which means that this instability is mediated by modes propagating at nearly the shear 
wave-speed in the direction of sliding. The propagation direction is determined by the positive 
sign in the last expression. It is important to stress in this context that 2 ; = 1 is not a solution of 


Eq. (34), again demonstrating the symmetry breaking induced by frictional sliding. 


The dilatational wave dominated instability exists for all physically relevant values of A, i.e. for 
0 < A < 1. Is it true also for the shear wave dominated instability? To address this question, we 


interpret A in Eq. (47) as a function of r = 7 /, parameterized by [3 and /. A(r) is a non-monotonic 


function which attains a maximum at 


^ 2(1 - /32) + /2 _ 2V(1-/32)(1-/I2 + /2) 

p 


(48) 


For realistic values of the friction coefficient / (i.e. /~0.2 —0.75), we have A*^™) which shows 
that the shear wave dominated instability is characterized by a small A. Furthermore, A(r) in 
Eq. (47) vanishes at T = Pl{l — p + P), which is also typically small due to the smallness of 
p. In fact, if we assume a small T and invoke a parabolic approximation for A(r) in Eq. (47), 


we obtain for the maximum A*^™) ~ /^/[4(1 — P)], which is just the leading contribution in the 
expansion of Eq. (48) in terms of p. We thus conclude that the shear wave dominated instability 


is localized in a relatively small region near the origin in the A —T plane. 

One implication of the above discussion is that since in the stability boundary of the dilatational 


wave dominated instability A increases linearly with T = 7 /, cf. Eq. (43), the shear and dilatational 
waves instabilities coexist only in a relatively small range of A’s, 0 < A < To explicitly 

demonstrate this, 3?(A) is plotted vs. q in Fig. for the two types of instabilities and various small 
A’s. We observe that indeed the two instabilities coexist for 0 < A < A^™'^, but only the dilatational 
one exists for A > A(”*) (see figure caption for details), and that the growth rate of the dilatational 
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Figure 4: K(A) vs. wave-number q for both the dilatational and shear wave dominated instabilities for various A’s, 
and 7 = 0.25. For the parameters used (together with / = 0.3 and /3 = 0.5), Eqs. (47l-(48l imply ~0.234 (marked 
with a vertical arrow) and 0.024. We observe that for A< A^™'^ the two types of instability coexist, while for 

A > A^™'^ only the dilatational one exists. Furthermore, the figure quantitatively verifies the prediction for the critical 
wave-number qi ”'’, demonstrating that indeed the shear wave dominated instability appears at a finite wave-number, 
unlike the dilatational wave one. The growth rate of the dilatational wave-like instability is larger than the one 
corresponding to the shear wave-like instability. 


instability is larger than that of the shear one. Furthermore, we see that indeed the dilatational 
wave dominated instability appears at a vanishing wave-number, while the shear wave dominated 
instability appears at a finite wave-number. The results for the shear wave dominated instability 
presented in Fig. were obtained numerically. We could have followed a similar procedure to 
the one taken in great detail in Sect. |6.1| and derive analytic results by systematically expanding 
around z = —l in the complex-plane. In order not to further complicate the presentation, we do 
not present this analysis here, but rather present numerical demonstrations of the main physical 
points. 

The analysis of the spectrum in the large kH limit presented so far has revealed two classes of 
elastodynamic-controlled unstable modes, one mediated by dilatational wave-like modes propagat¬ 
ing in the direction opposite to the sliding motion (corresponding to solutions near z = l//3) and 
one mediated by shear wave-like modes propagating in the direction of sliding (corresponding to 
solutions near z = —l). Related observations on the directionality of unstable modes and rupture 


along bi-material frictional interfaces have been previously made, see for example Ranjith and 


Ricej (i 

1001); Cochard and Rice ( 

2000); Adams (2 

000); Weertman (1980); 

Andrews and Ben-Zion 

(1997 

Ben-Zion and Huang (20C 

2); Adams (1995 

1998); Harris and Day 

(1997); Xiaet al. (200% 


Ampuero and Ben-Zion (2008). 


Finally, we emphasize that one cannot naturally superimpose the results of Eq. ([^ and Fig. [2] 
on those appearing in Fig. in a generic manner because the former results depend on while 
the latter do not (i.e. they are valid for H~^). In particular, the growth rate in Fig. [^decays 
as and its relative magnitude compared to the growth rates in Fig. [^depends on the value 
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of H. It is important to note, though, that for a real system with a given H our results allow one 
to calculate the growth rate of all of the instabilities discussed above and determine which is the 
largest. 


6.3. The quasi-static limit 

Up to now we have found two classes of elastodynamic-controlled instabilities, one related to 
dilatational waves and one to shear waves. In addition to these, there exists also a quasi-static class 
of unstable modes at extremely small A’s, which is qualitatively different as it is not elastodynamic 


in nature. This quasi-static instability has been discussed quite extensively in Rice et al. (2001) 


where the analysis has been performed for general bi-material interfaces (i.e. not only for a strong 
contrast). Our goal here is just to briefly summarize those results of Rice et al. (2001) which are 
relevant to our discussion. In order to see how the quasi-static limit emerges, we need to take the 
limit of large wave-speeds Cg^d —)■ oo in the linear stability spectrum in Eq. (20), while keeping their 
ratio 13 fixecQ Obviously, the friction part remains unaffected and the elasticity parts are changed 
according to 

1 17 1 A(l — u) 2(1 — 2u) 

kd ^ \k\, 9l ^ -JTT; 92 —^ - 77 — • (49) 


3 - 4.12 


3-4i2 


The resulting quasi-static linear stability spectrum can be analyzed following Rice et al. 
leading to the stability condition 


( 2001 ) 


(50) 


where an instability appears when the left-hand-side is smaller than the right-hand-sid^ We 
do not report here on the detailed calculations that show that the leading elastodynamic effect 
on this instability branch enters only to quadratic order in T = 7 / (which, as discussed above, 
quantifies the importance of elastodynamic effects) and we do not explore the range of existence 
of this instability branch with increasing T. For our purposes here it would be sufficient to note 
that to leading order, the instability condition in Eq. (50) is independent of T. The quasi-static 


instability emerges at a finite wave-number (again, this is not shown explicitly here), >0 and 
consequently, a finite range of unstable modes exists above the threshold. 

The quasi-static instability exists at extremely small values of A, typically of the order of 
10“^ due to the appearance of higher powers of / and j3 (both smaller than unity) in Eq. (50). 
Yet, since the stability condition for both the dilatational and shear wave dominated instabilities 
satisfies AocT for sufficiently small T, there exists a small region near the origin of the A—T plane, 
where only a quasi-static instability can be found. All of these issues will be addressed next, where 
we construct the stability phase diagram of the problem in the large kH limit. 

6 . 4 . Stability phase diagram 

One of the hallmarks of a linear stability analysis is a stability phase diagram in the space of 
the relevant physical parameters. The detailed analysis presented above allows us at this point to 


®Note that if the limit Cs,d^oo is taken at the level of the fields themselves in Eq. (131, the solutions become 

degenerate and additional independent solutions should be included. _ 

Equation 


( |50| coincides with the last equation on page 1890 of Rice et al. (20011 (the equations in that paper 
are not numbered). Note, however, that f} in Rice et al. ( 2001| denotes the Dundurs parameter (which is a function 
of the 4 linear isotropic elastic moduli of the two sliding bodies) and is different from our /3. In the limit of infinite 
material contrast, which is the limit we consider, the Dundurs parameter simply equals \/~3 (where /3 is defined in 
Eq. (23l). One should bear this in mind when comparing Eq. (|50[) to the results of Rice et al. (20011. 
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analytically construct such a stability phase diagram. We have identified three classes of insta¬ 
bilities in the large kH limit (a dilatational wave dominated instability, a shear wave dominated 
instability and a quasi-static instability) and the corresponding stability boundaries are given in 
Eqs. (|4^, (|47|) and ([50|). We rewrite these as 


Arf(r) = r /3 (1//32 - 2), A,(r) = r - 


r2(i -,02^ 
^(i-r)/2 ’ 


A,,(r)~i/2/34 + o(r2) , (51) 


where we added the subscripts d, s and qs to correspond to “dilatational”, “shear” and “quasi¬ 
static” instabilities, respectively. 

We treat all of these (in)stability boundaries as functions of F = 7 /, parameterized by / and {3. 
While there is some degree of arbitrariness in this choice, we strongly believe that it is the most 
natural way to represent the interplay between the various physical effects in the problem, within 
the framework of a two-dimensional phase diagram. An instability is implied whenever A is below 
at least one of the Ac;(r), As(r) or Ag 5 (r) lines in the A —F plane. Note that while the results 
for Arf and A^ are exact, the one for Aqg is given to leading order in F. This will be enough for 
our purposes here. Finally, we stress that the phase diagram to be presented and discussed below 
pertains to the large kH limit; in the limit of small kH, as discussed in Sect. homogeneous 
sliding is unconditionally unstable. 

As there exist various types of instabilities, one is interested in understanding under what 
conditions they coexist and what is the upper stability boundary, i.e. the line which separates the 
region in parameter space where no instabilities exist at all from the region in which at least one 
instability exists. In Fig. the A —F stability phase diagram for two sets of values of (/,/3) is 
shown. Note that negative values of A, which correspond to steady-state velocity-weakening, are 
not shown as in this case there is an instability independently of other parameters. 

Let us discuss in detail the linear stability phase diagram. A(i(F) is a straight line spanning 
the whole range of A values, 0<A<1, while the other lines are localized near the origin (readers 
are advised to remind themselves the physical meaning of the dimensionless parameters used here, 
as discussed around Fq. (|21[)). In this sense, the dilatational wave dominated instability is the 
main instability mode of the system in the large kH limit. As discussed in Sect. |6.1[ the trends 
in A£;(F) are clear; for a fixed F, increasing A (essentially increasing d^f, recall that this stability 


boundary is independent of the direct effect dvf, as discussed around Fq. (44)) promotes stability. 
Alternatively, for a fixed A, instability is promoted by increasing F (either by enhancing the 
elastodynamic bi-material effect ~^/cs or by decreasing the normal load ao). 

As(F), as discussed in Sect. |6.2[ is a non-monotonic function that is bounded in the region 
0 < A < A^™'^ and crosses zero at a finite small F, F = /2/(l — + /^) (recall that A^™) <C 1 


is defined in Fq. (48)). This implies that, since Ac;(F) is a straight line starting at the origin, 
the shear and dilatational wave dominated instabilities always have a range of coexistence. The 
question then is whether there exists a parameter range where the shear wave instability exists 
and the dilatational one does not. It is evident from Fq. © that this depends on the value of 
j3. For small F, we have As(F)~F -|- O (F2), i.e. As(F) starts linearly with a unit slope. Ad(F) is 
always linear with a slope /3 (I //32 — 2). Therefore, for /3 — 2) >1 the shear wave dominated 

instability coexist with the dilatational one, while for fi (l//3^ — 2) < 1, there exists a region above 
the Arf(F) line in which the shear wave dominated instability exists, but the dilatational one does 
not. In the upper panels of Fig. [^we used /3 = 0.45, which is an example of the former, while in 
the lower panels we used /3 = 0.6, which is an example of the latter. 
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Figure 5: Stability phase diagram in the A—F plane. Each row presents the same data (corresponding to the same / 
and P), but the left column is in linear coordinates and the right one in logarithmic coordinates. Sliding is unstable 
beneath the upper stability boundary (dashed black line), which is the upper envelope of the three stability lines of 
Eq. (511. The shaded colored areas correspond to regions of the phase diagram in which one or more instabilities 
exist (i.e. below the upper stability boundary, where each color corresponds to a different instability). In the two 
right panels, all dilatational modes are unstable in the region to the right of the dotted line, see Eq. (451. The 
quasi-static instability exists in a very small region near the origin of the A —T plane and hence is seen only in the 
logarithmic plots (right panels). Note also that since we only consider the quasi-static stability boundary to leading 
order in T, we truncated it at small values of T. For concreteness, we chose to truncate it when it intersects one of 
the other stability boundaries, though there is nothing physically special about these points (beyond the fact that 
they belong to the upper stability limit). 


Aqs(r) is, to leading order, a positive constant independent of F. Since both Arf(r) and As(r) 
vanish linearly at small F, there always exists an instability region controlled by quasi-static modes. 
This quasi-static stability boundary joins either the Arf(F) line or the As(F) line, depending on the 
value of fj. It happens at F~/^/3^/[4(l//3^ — 2)] in the former case and at F~/^/3^/4 in the latter. 
The two possibilities are shown on the two right panels in Fig. The upper stability boundary 
always starts with Aq 5 (F), and then either merges directly with A(i(F) (if /3 (l//3^ — 2) >1) or 
first merges with As(F) (if /3 (l//3^ — 2) < 1) which then merges with A(i(F). One way or the 
other, except for a small region near the origin of the A—F plane, the upper stability boundary is 
determined by A^(F). As mentioned above, we do not consider higher order corrections to Ags(F) 
in terms of F and hence we truncate the latter when it intersects either As(F) or A(i(F) in Fig. 
(we stress, though, that there is nothing special about the intersection point from the perspective 
of Aqs(F)). Finally, note that in a region of coexistence, the instability with the largest growth 
rate K(A) will be observed physically, cf. Fig. 
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With this we complete the discussion of the linear stability analysis in the large kH limit for 
large contrast bi-material interfaces described by standard rate-and-state friction. In what follows, 
we will discussed some generalized rate-and-state friction models, which incorporate a modified 
non-instantaneous response to normal stress variations. 


7. Generalized friction models: Modified response to normal stress variations 


Up to now we considered the standard rate-and-state friction model in which the frictional 
resistance does not exhibit a finite time response to variations in the local normal stress. That 
is, the friction stress T(x,t) was affected by the local interfacial normal stress cr{x,t) only through 


Eq. ([TT|), and the response was instantaneous. Nevertheless, some experimental work ( 

Linker 

and Dieterich 1992| Dieterich and Linker, 1992t Prakash and Clifton 1992[ 1993; Prakash 

1998t 

Richardson and Marone 


Bureau et al. 

2000 

) indicated that the frictional resistance may 


exhibit a finite time response to normal stress variations. To take this possibility into account 
we adopt here two experimentally-based modifications of the constitutive relation and incorporate 
them into a generalized analysis. The goal is to understand the physical effects of the modified 
constitutive relations on the stability analysis. 


The first modification we consider is due to Prakash and Clifton (1992, 1993); Prakash (1998) 


and amounts to taking the time scale T in Eq. (10) to be finite. A hnite T means that when 


the interfacial normal stress cr(x, t) varies, the friction stress t{x, t) does not immediately follow 
it as in Eq. 0- T is commonly expressed as T = ap^D/v, where is a dimensionless and a 
positive parameter which measures T in units of the already existing time scale in the model, D/v. 
Consequently, F{T,a,v,(j)) in Eq. ([^ takes the form 


f = F(t, (t,u,(/>) = - 




(r - cr/((/>,u)) . 


(52) 


The instantaneous response in Eq. 0 is recovered in the limit Op^ 0. Eor any ap^ > 0, 
the response is gradual such that the larger ap^^ the slower the response. As the variation of 
the interfacial normal stress cr{x, t) with slip (the elastodynamic bi-material effect) is the major 
destabilizing effect in the analysis presented up to now and since ap^ > 0 delays the destabiliz¬ 
ing effect on the friction stress T{x,t), we expect Op^ > 0 to promote stability. The finite-time 


response to normal stress variations in Eq. (52) has been invoked by several authors in relation 


to the regularization of the Coulomb friction law, for example in the context of the stability of 


steady frictional sliding (Ranjith and Rice, 2001), which was already mentioned above, and dy¬ 


namic rupture propagation (Cochard and Rice, 2000; Ben-Zion and Huang, 2002) along bi-material 
interfaces. 


The second modification we consider is due to Linker and Dieterich (Linker and Dieterich 


1992, Dieterich and Linker, 1992), who interpreted their experiments as suggesting that the state 


variable (j) is directly affected by variations in the normal stress. In particular, they proposed that 
G(r, a,v,4>) in Eq. (1^ takes the form 


= —-a„— , 


(53) 


where is a positive dimensionless parameter. When =0, Eq. is recovered with g{x) = 
1 — X (note that ^'(1) = —1). To qualitatively understand the effect of > 0 on the frictional 
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stability, we consider a fast variation in the normal stress (i.e. a large Idl) such that the last term 


on the right-hand-side of Eq. (53) dominates the other two terms. In this case, we can eliminate 
the time derivative to obtain 


(5(/> ~ 


crod^f 


o-Q^/ - -arD^cr , 


(54) 


where we set a = ao. 

The last result may appear somewhat counterintuitive. To better understand it, note that for 
apc =0 (which is assumed here to simplify things), (jQ6f contributes to the variation of the friction 
stress St as in Eq. (12). According to Eq. (|54|), when the normal stress reduces, ci < 0, we have 


aoSf>0 which means the latter makes a positive contribution to the friction stress (for > 0 ). 
This appears to suggest that the bi-material effect in this case is stabilizing. This is not quite the 
case, because the effect of 5a on St contains in fact also the last term on the right-hand-side of 
Eq. (12). Together, we obtain dr ~ (/ — Q;^^)dcj (recall that Sayy = —Sa), which shows that as 


long as / > Opp,, a reduction in the normal stress still leads to a reduction in dr, i.e. it remains 
a destabilizing effect (later on we will see that this is not a strict stability condition). Yet the 
magnitude of the destabilizing effect is reduced when > 0 (i.e. it is determined by / — 


instead of / alone), which indicates that the modihcation in Eq. (53), with > 0, promotes 
stability. Consequently, the simple considerations discussed here suggest that replacing Eqs. 
and 0 with Eqs. ([52|)-([5^, for Op^ >0, will facilitate stability. 

Next, we aim at rigorously studying the generalized model, which incorporates the modified 


response to normal stress variations in Eqs. (52)-(53), in the large kH limit. For that aim, we first 


derive the linear stability spectrum for this case. Contrary to the analysis in Sect. we should 
explicitly include now perturbations of the friction stress r(x, t) =ro-|-A,- exp[At — ikx] since r(x, t) 
satisfies a dynamical equation of its own. The linear stability spectrum reads 


f^kdgi{l + Opp^A— ) -igkfg2 + 

Vo 


ik pL Opp, ^2 A (To A 


A+^ 




A+^ 


(Adpf + ^dyf^ = 0 , (55) 


which reduces to Eq. (20) for = app, = 0, |fl'^(l)| = 1 and gi ^2 —^ Gi, 2 - The non-dimensional 
spectrum takes the form 


7(1 - iqz) gi{z,l3) (1 -iapp.qz)- ifg2{z,(3) A'yqUpp g2{z,l3)z- iz {A - iqz) = 0 , 

(56) 

which reduces to Eq. (34) when = 0. Our next goal would be to analyze this linear 

stability spectrum. 

7.1. Dilatational wave dominated instability 


To analyze Eq. (56), we closely follow the procedure described in Sect. 6.1 We first look for 


solutions located near z = l/fd, corresponding to a dilatational wave dominated instability. The 


expansions in Eqs. (36) and (^39^ remain valid and can be substituted into Eq. (56) to yield 


Kd ^ 


if'yjl -iq/l3) (2- l//3^) + i {A - iq//3)//3 - j qUpp (2- l//3^)/^ 
7 ( 1 - iq/P) (1 ^ /\/l//3^ - 1 - i Opc Q/P) //3^ T * 7 9 “ld \/1//3^ - l//3^ 


(57) 
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Here, as before, the stable solution corresponds to the minus sign and unstable one to the plus 
sign. Equation (57) is then substituted in Eq. (36) to obtain z, from which the growth rate is 
derived according to 3^(A) = ( 79 (z). 

While the resulting analytic expression for 3f?(A) is readily available, it is a bit lengthy and we do 
not report it explicitly here. Yet, some analytic insight can be gained directly from Eq. (57). Eirst, 
we note that always appears through the combination / — ; in the numerator it appears 

through 'jq{f — a^^){2 — l//3^)//3 and in the denominator through i'jqif — a^j^)jl3^. 
This is in line with the qualitative discussion above, which indicated that the main effect of is 
to reduce the effective friction coefficient. However, it is crucial to understand that while enters 
the problem through the combination / — , / also appears independently (cf. the first term in 

the numerator of Eq. (|57|)). Furthermore, while / — is always multiplied by the wave-number 
q, f appears independently of q. This structure will have direct implications for the stability 
boundary, which corresponds to g—5-0, as will be discussed below. 

The other new parameter, Op^^, is also multiplied by q in Eq. (57). In fact, it introduces a new 
term proportional to oip^q^ in the denominator of Eq. (57), which does not exist in the theory 
with Op^ = 0. This has interesting implications for the behavior of the growth rate as q^ oo. 
Our previous analysis for Op^ = 0 showed that 3^(A) approaches a finite constant as q^ oo. For 
Qfpp > 0, the presence of the new term proportional to otpQq^ in the denominator of Eq. (57) and 
the fact that the largest power of q in the numerator is linear, implies that 5?(A) —?■ 0 as g— t-oo. 
This suggests that even if all wave-numbers are unstable, the growth rate vanishes for sufficiently 
large q. This property of the generalized model is appealing from a basic physics perspective and 
as such it constitutes an improvement relative to the standard rate-and-state friction model. 

The analytic properties discussed above are demonstrated in Fig. where we show 3?(A)((?) as 
obtained from a direct numerical solution of Eq. (56) and the analytic prediction (3^(A) = g9'(z), 
where z is given in Eq. (36) with Kd of Eq. (57)) for two sets of the parameters (A, /, 13, 7 ), =0 

and various values of Up^ (only 3f?(A) > 0 is shown). That is, we isolate the effect of Op^ on the 
stability problem. A few conclusions can be drawn from the figure. First, we observe that the 
analytic prediction is in a very good quantitative agreement with the exact numerical solution for 
all of the parameters considered, yet again lending strong support to the theoretical approach. 
Second, as expected, we observe that increasing > 0 reduces the range of instability and the 
magnitude of the growth rate (i.e. it promotes stability). Finally, we observe (right panel) that 
as predicted analytically, I?(A) decays to zero at large g’s even when all q's are unstable. We also 
observe that an infinite range of unstable modes can become hnite upon increasing otp^. 


In Fig. 
as in Fig. 


we elucidate the effect of (X „ on the stability problem. We follow the same scheme 


but now set = 0 and vary a^p. We again observe that the analytic prediction 
provides a good quantitative approximation to the exact numerical solution (some deviations are 
observed with increasing a^p in the left panel). We also observe that as the simplified analysis 
above predicted, increasing a^p promotes stability. Note that since we consider here Op^ =0, 3^(A) 
approaches a constant as g —)• 00 because a finite a, „ does not introduce higher order powers of q 


into the expression for Kd in Eq. (57), as compared to the a^p =0 case discussed in Sect. 6.1 As 


Fig. [^clearly shows (right panel), the value and sign of the constant does depend on a^p. Finally, 
we note that an instability exists also when / — a^p < 0 (see left panel), demonstrating that a^p > f 
does not eliminate altogether the destabilizing elastodynamic bi-material effect, as implied by the 
simplified analysis at the beginning of this section. 

The last issue we need to discuss is the threshold condition (stability boundary) associated with 
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Figure 6: The dimensionless growth rate 5R(A) vs. wave-number q from a direct numerical solution of Eq. (561 (thick 
lines) and the analytic prediction corresponding to Eq. ( |36| l with K,d of Eq. (571 (thin lines). Here we set =0 
and use various values of Opp. In the left panel we use 7 = 3.2 and in the right one 7 = 4 such that E = 0.96 and 1.2, 
respectively (the rest of the parameters are as in Fig.j^ i.e. / = 0.3, /3 = 0.5 and A = 0.7.). 




Figure 7: The same as Fig. but with =0 and varying values of a^p,. 


the dilatational wave dominated instability. As we have seen in relation to Eq. ( |40| ), the threshold 
condition emerges from the sign of the g-independent part of the numerator of '^[nd] in the latter 
equation. On the other hand, as discussed above, Op^ and in Eq. (57) are multiplied by q, 
which implies that they do not affect the threshold condition. Consequently, Eq. (43) remains 
the threshold condition for the dilatational wave dominated instability in the generalized model, 
independently of and Opp. In the next subsection, we discuss the threshold condition associated 
with the shear wave dominated instability within the generalized model. 


1.2. Shear wave dominated instability 

We now consider the shear wave dominated instability, 
z = — l in the complex plane. Repeating the analysis that led 
and the stability threshold read 


which corresponds to solutions near 
to Eq. (47), the critical wave-number 


1 

L 

l-7(/-app) 

l-7(/-app) 

a 

LXpc 

_2app7^1 -/32 

2apc7\/l - 


and A = 7 /- 7 v^W^(l-bQ^p)g(-) . 


(58) 
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We first observe that the stability threshold, unlike the one for the dilatational wave dominated 
instability, depends explicitly on both a, 


ipp and (the latter only through This happens 


because the instability occurs at a finite wave-number. Equation (58) shows that while depends 
on / and only through the combination / — (in addition to its dependence on Op^), 
threshold condition depends both on / and / — . 


the 


The limit Op^ —)• 0 should be handled with some care as a few terms in Eq. (58) appear to diverge 
upon a naive substitution of =0. Nevertheless, this limit exists whenever 1 — 7 (/ — >0 

and takes the form 


iM = 


13 ^ 

1 -7(/ - Opp) 


and 


A = 7/ - 


7 ^( 1 -/ 3 ^) 

1 -7(/ - app) 


(59) 


Einally, Eq. (47) is recovered once (Xp = 0 is substituted in the above expressions. With this we 


complete the large kH analysis of the generalized models in the elastodynamic regime (we do not 
consider here the quasi-static limit). 


8. Brief discussion and concluding remarks 

In this paper we considered in quite general terms the linear stability of homogeneous sliding 
along frictional interfaces separating strongly dissimilar elastic materials, with a focus on finite size 
effects, elastodynamic effects and velocity-strengthening friction. The linear stability problem is 
studied, analytically for the most part, within the constitutive framework of generalized rate-and- 
state friction models, including velocity-weakening associated with the maturity/age of contact 
asperities, instantaneous rheological strengthening, steady-state velocity-strengthening and a reg¬ 
ularized response to normal stress variations. We considered finite size systems, of height H, and 
analyzed the stability spectrum in both the small and large kH limits. The various competing phys¬ 
ical effects, most notably the destabilizing elastodynamic bi-material effect and stabilizing effects 
associated with the interfacial constitutive behavior, are quantified through several dimensionless 
parameters. 

We showed that there exists a universal instability (independent of the details of the friction law, 
but possibly dependent on the assumed large material contrast) characterized by a wave-number 
k^H~^ and a maximal growth rate 51?[A] oc/cs/77, mediated by waveguide-like modes. The role of 
boundary conditions for finite size systems has been also highlighted. Eor large systems, H ^oo, we 
provided a comprehensive quantitative picture of the stability phase diagram. We showed that in 


addition to the previously derived quasi-static instability modes (Rice et al., 2001), which exist at 


sufficiently small sliding velocities, there exist also dilatational wave dominated instability modes 
propagating in the opposite direction to the sliding direction and shear wave dominated instability 
modes propagating in the sliding direction. The former appear to be the dominant instability 
modes over a broad range of physical parameters. Our results allow to determine, for every set 
of physical parameters relevant to a specific frictional system, which of the discussed instabilities 
features the largest growth rate. In a certain parameter range the instability is manifested through 
unstable modes at all wave-numbers, yet the frictional response is shown to be mathematically 
well-posed. The stabilizing roles played by a regularized response to normal stress variations are 
quantitatively accounted for. 

All in all, our analysis shows that steady sliding along strong bi-material frictional interfaces 
is quite generically unstable, even in the presence of steady-state velocity-strengthening friction of 
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a general form. The corresponding problem for steady-state velocity-weakening friction is known 
to be generically unstable, but this instability is of different origin compared to the instabilities 
discussed in this paper. The discussed instabilities depend crucially on the destabilizing bi-material 
effect and consequently these generic instabilities are not expected to persist in frictional interfaces 
separating identical materials. Our results should be relevant to any frictional system exhibiting 
strong material contrast along the frictional interface (and possibly also to finite contrast inter¬ 
faces) and steady-state velocity-strengthening, for example geophysical systems such as mature 
earthquake faults ( Ben-Zion| 2001, 2008) and engineering/tribological systems such as an elastic 
brake pad sliding on a rigid substrate (jBehrendt et ah' 


MITj ). 


The results of the linear stability analysis are possibly related to rupture dynamics along bi¬ 


material interfaces (Ben-Zion 2001,2008). While the exact relations can be elucidated by following 


the instabilities into the nonlinear regime, most probably only numerically, one can speculate about 
these relations. For example, the dilatational wave dominated instability, which propagates in 
the direction opposite to the sliding motion and emerges at A: —)• 0 (i.e. large wavelength), may 
correspond to crack-like (non-localized) super-shear rupture fronts propagating in the so-called 
“un-preferred direction” (i.e. opposite to the slip direction of the more compliant material, for 
any material contrast). Likewise, the shear wave dominated instability, which propagates in the 
direction of sliding motion and emerges at a hnite k, may correspond to pulse-like (localized) 
shear rupture fronts propagating in the so-called “preferred direction” (i.e. the slip direction of 
the more compliant material, for any material contrast). These possible relations can be tested 
exp er iment ally. 

In the future, it would be interesting to quantify how the existence of a finite bi-material contrast 
affects the stability analysis presented here for the large contrast limit, where the destabilizing 
elastodynamic bi-material effect is the largest. Furthermore, other theoretical issues should be 
considered; for example, as the dilatational and shear wave dominated instabilities propagate with 


high velocities in the lab frame of reference, issues of absolute vs. convective instabilities (Lifshitz 


and Pitaevskii 1981; Huerre and Monkewitz, 1990) may become important. Moreover, in the 


presence of lateral boundaries (i.e. for realistic systems of finite size L in the x-direction) the 


so-called global instability (Lifshitz and Pitaevskii, 1981; Huerre and Monkewitz, 1990) may be 


relevant due to the lateral boundary conditions which induce a coupling between different linear 
propagating modes. Finally, the finite size L naturally implies a modification of some of the small 
wave-number instabilities discussed above. 
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Appendix A. Derivation of Eq. (30) 


We start by expanding first the basic functions appearing in Eq. (24) (see also Eq. ©) as 

(Appendix A.l) 


kd{A,k) ~ ±— , 
Cd 

,,,,,, , H5K 

coth(iffcrf) ~ ±- 


fes(A, k) ~ ±— , 

Cc 


Cd 


tanh{Hks) — ± tanh(ffAo/cs 


where the ± correspond to the stable/unstable branches, respectively. As will be seen right away, 
the distinction between the two branches disappears in the final expression. Using these expansions, 
we obtain the following relevant combinations 


ks(A,k) tanh.(H ks) — tanh(ffAo/cs) , 

Cs 

ks{A, k)kd{A, k) ~ , 

coth{H kd) taiiih{H ks) ^ tanh(ffAo/cs) . 

Cd 


(Appendix A.2) 


These expressions allow us to expand the three contributions to S{A,k,H) (after dropping the 
common prefactor /i). 

First, we consider the contribution to S{A,k,H) in Eq. (24), associated with Sayy, 

GiiA,k,Hy 


ik f G'2(A, k, H) = —ik f [2 — 


coth.{Hkd) 


ik f 


AVc 


= -ikf\^2- 
ifcsk 




kdks coth(Fffcd) tanh{Hks) — k'^ 
fcsk 


kdks coth.{Hkd) tan]i{Hks) tanh(FfAo/cs)(5A 


= ± 


(3‘^H tan(i7| Ao|/cs)(5A 

(Appendix A.3) 


It is important to note that the latter depends on the ratio of the two small quantities k and 6A 
and hence contributes to zeroth order. Therefore, in what follows we can expand quantities to 
zeroth order. 

For the contribution associated with 5axy, we have 


kd{A,k)GiiA,k,H) = 


kd kd) A"^/c^ 


AVc 


|Ao| 


kdks cotli{Hkd) tanh{Hks) — k"^ ks ta'nh.{Hks 


Cs tan(Ff|Ao|/cs) 

(Appendix A.4) 


where in the last step we used the fact that Aq is purely imaginary, which implies that tanh(FfAo/c^) = 
±itan(i?|Ao|/cs) (± here correspond to Aq = ±z|Ao|). We thus conclude that the contribution 
corresponding to 6axy is finite and independent of the sign of A(Ao). Note that we started by 
expanding around (A = AQ,k = 0), which corresponds to 6axy{A,k —)■ 0) = 0 (see Eq. (25) and 
the discussion around it), and now we find that 6axy{A —)• AQ,k —)■ 0) is actually finite. There 
is, however, no contradiction here. The point is that previously we used 5axy{A,k) = 0 to derive 
the dispersion relation Awg{k) and then took the A: —)• 0 limit, and here we treat 5axy{A,k) as a 
function of two independent variables, where A(A;) is not known (the goal is to determine it from 
S{A,k,H) = 0). 
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Finally, the friction law contribution is simply evaluated at A = Ao, i.e. 

A A(A) + A Aq Aq + a |Ao| / —|Ao|(1 — A) ± i(A + |AoP) 


Cs a l) 7(^0 A 1) 


l + |Ao|2 


I (Appendix A.5) 


where Aq = arid again the ± correspond to Ao = ±i|Ao|. Collecting all three contributions 

we end up with Eq. (30) of the main text. 


Appendix B. Finite height analysis with imposed shear stress 

The analysis of finite size systems (i.e. in small kH limit) in Sect, [^was performed under an 
imposed tangential velocity boundary condition at y = H. As mentioned above, we expect the type 
of boundary conditions at y = H to affect the stability problem. To demonstrate this, we consider 
two situations which are identically the same (in terms of both parameters and geometry) except 
that in one case we impose a tangential velocity vq at y = H (i.e. as on the left of Eq. ([^, which 
was adopted throughout the manuscript until now) and in the other case we impose a shear stress 
tq at y = H (i.e. as on the right of Eq. ([^). Here vq and tq are related by the friction law under 
steady-state conditions. The spectrum equation for the imposed shear stress case is derived along 
the same lines as for the impose tangential velocity, but the result is a bit too lengthy to be reported 
on explicitly here. Furthermore, we restrict ourselves in this Appendix to its numerical analysis, 
though in principle analytic progress in the spirit of the analysis performed in the manuscript 


can be pursued. In Fig. B.8 we plot the normalized growth rate H^[A]/cs as a function of the 
normalized wave-number kH for the two cases. It is evident that the result quantitatively depends 
on the type of boundary condition, though no marked qualitative differences are observed. 



Figure B.8: The growth rate 5R[A] (in units of Cs/H) vs. Hk, for small Hk, as in the left panel of Fig. The 
upper (red) curve corresponds to velocity-controlled boundary conditions (it already appeared in the left panel of 
Fig. HI with jf— 0.9 and n = 0, though the range of Hk is smaller here). The lower (green) curve corresponds to 
stress-controlled boundary conditions (see text for details). 


Appendix C. Discontinuities (gaps) in the spectrum: Non-localized modes 

The full numerical solution of the implicit linear stability spectrum in the kH^oo limit, shown 
on Fig. exhibits a discontinuity (gap) at the transition from the unstable to the stable part of 
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the solution. Our goal here is to propose an explanation for the origin of the gap in the spectrum. 
As a background, note that while in Sect. the H^oo limit was taken directly at the level of the 
linear stability spectrum of Eq. (24), cf. Eq. (p^ , one could alternatively take the limit at the level 


of the elastodynamic solutions of Eqs. (13). In this case, we demand that the perturbations decay 
away from the interface, i.e. as y—^oo, which implies that A3 = A4 = 0 and only Ai ^2 remain finite. 
Imposing the two boundary conditions at the interface, y = 0, yields the instability spectrum in 


Eq. (34). 


In this H ^ oo limit, only the decaying solutions {exp[—exp[—are physically rel¬ 
evant. We will now show that the gap in the spectrum is related to an increasing solution, 
which is non-physical in the strict H ^ oo limit. In particular, we consider a solution involv¬ 
ing {exp[ksy],exp[-kdy]}, that is 


6ux{x,y,t) 

6uy{x,y,t) 


k 

-ikd 


kg 

ik 


.42exp|-M\ [ ] 

Asexpp.y] ) ^ 


(Appendix C.l) 


By imposing the boundary conditions at y = 0 and properly nondimensionalizing all relevant 
quantities we obtain the following linear stability spectrum equation 


s{z, q) = ^ — iqz) —■sj\ — gi{z, /3) — ifg 2 {z, (3) — (A — zyz) =((Appendix C.2) 


where 


gi{z,l3) = - 


1 -|- Vl — — 13‘^Z^ 


and 


g2{z,l3) = 2 +gi{z, ^^.^ppendix C.3) 


which are the counterparts of Eqs. (34)-(35). 



Figure C.9: The dimensionless growth rate 5R(A)=qQ'(«) vs. the dimensionless wave-number q, as in the left panel 
of Fig.j^ The solid (red) line is identical to the 7 /= 0.96 curve in the left panel of Fig. The dashed (green) curve 
corresponds to solution of the spectrum in Eq. ([Appendix C.2[| for the same set of parameters. 


In Eig. C.9 we present numerical solutions of Eqs. (34) and (Appendix C.2) for the very same 
set of parameters. It is observed that while each of these solutions is discontinuous in itself, 
superposing the two solutions generates two continuous functions, each of which is composed of 
two segments from different solutions. We thus propose that the gap observed in Eig. is related 
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to the solution in Eq. ([Appendix C.l|), 

y- 


which features a non-physical exponential divergence as 


•oo. In a real system, with a finite height H (however large), increasing solutions always exist 
and will naturally lead to the regularization of the discontinuities that emerge in the strict H^oo 
limit. 
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